
##TMESCORE
#引用包
library(limma)
library(estimate)
inputFile="OV_TPM.txt"       #表达输入文件
setwd("E:\\1.Ovarian_Treg\\1.picture\\14.TME")       #设置工作目录

#读取文件,并对输入文件进行整理
rt=read.table(inputFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)

#删除正常，只保留肿瘤样品
group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
data=data[,group==0]
out=data[rowMeans(data)>0,]
out=rbind(ID=colnames(out),out)

#输出整理后的矩阵文件
write.table(out,file="uniq.symbol.txt",sep="\t",quote=F,col.names=F)

#运行estimate包
filterCommonGenes(input.f="uniq.symbol.txt", 
                  output.f="commonGenes.gct", 
                  id="GeneSymbol")

estimateScore(input.ds = "commonGenes.gct",
              output.ds="estimateScore.gct")

#输出每个样品的打分
scores=read.table("estimateScore.gct", skip=2, header=T)
rownames(scores)=scores[,1]
scores=scores[,3:ncol(scores)]
colnames(scores)=gsub("\\.", "\\-", colnames(scores))
colnames(scores)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*", 
                      "\\1\\-\\2\\-\\3", colnames(scores)) 
scores=scores[,!duplicated(colnames(scores))]
scores=t(scores)
out=rbind(ID=colnames(scores), scores)
write.table(out, file="TMEscores.txt", 
            sep="\t", quote=F, col.names=F)

#TME
###相关性散点图
library(limma)
library(estimate)
library(Hmisc)
library(tidyverse) 
risk=read.table("rs_tcga.txt",header = T,
                row.names = 1,sep = '\t',check.names = F)

sameSample=intersect(row.names(risk), row.names(scores))
risk=risk[sameSample,"riskScore",drop=F]
scores=scores[sameSample,,drop=F]
risk_TME=cbind(scores,risk)

save(risk_TME,file = "risk_TME.rdata")
data=risk_TME
##删除偏离值的样本
p_list <- list()
max(data$riskScore)
#data$riskScore=log2(data$riskScore)
for (i in 1:(ncol(data)-1)) {
  res <- rcorr(data$riskScore, data[,i])
  p_value <- signif(res$P[1,2], 2)
  cor_value <- round(res$r[1,2], 2)
  
  # 每次新建一个绘图数据框：
  data_new <- data[,c(ncol(data),i)]
  colnames(data_new) <- c("riskScore", "y")
  
  p <- ggplot(data_new,aes(x = riskScore, y = y))+
    geom_point(color = "#988d7b")+
    geom_smooth(method = "lm", formula = y ~ x, 
                # 调整置信区间颜色：
                fill = "#273359", color = "#00aeef", alpha = 0.8)+
    theme_bw()+
    ylab(colnames(data)[i])+
    theme(
      # 去除网格线：
      panel.grid = element_blank(),
      # 修改坐标轴标签
      axis.title = element_text(face = "bold.italic"),
      # 标题居中：
      plot.title = element_text(hjust = 0.5, size = 10)
    )+
    labs(title = paste0("r =", cor_value, ", q = ", p_value))
  p_list[[i]] <- p
}

# 拼图
library(cowplot)
library(patchwork)

p <- plot_grid(p_list[[1]], p_list[[2]], p_list[[3]],
               p_list[[4]],ncol = 2)

ggsave("risk_TME.pdf", plot = p, height = 5.5, width = 6)
##

#boxplot
#"StromalScore"  "ImmuneScore"   "ESTIMATEScore" "TumorPurity"  
#引用包
library(limma)
library(ggplot2)
library(ggpubr)
expFile="rs_tcga.txt"       #表达数据文件
gene="ImmuneScore"

#读取基因表达文件,并对数据进行处理
rt=read.table(expFile, header=T, sep="\t", check.names=F,row.names = 1)
rt$risk=ifelse(rt$riskScore>median(rt$riskScore),"high","low")
rt=rt[,"risk",drop=F]
risk_TME=risk_TME[,1:4]
data=merge(risk_TME,rt,by=0)
row.names(data)=data$Row.names
data=data[,-1]
#colnames(risk_TME)
data=data[,c(gene,"risk"),drop=F]
#临床相关性分析，输出图形结果

colnames(data)=c("gene", "clinical")
#data$gene=log2(data$gene+1)
#data=data[(data[,"clinical"]!="unknow"),]
#设置比较组
group=levels(factor(data$clinical))
data$clinical=factor(data$clinical, levels=group)
comp=combn(group,2)
my_comparisons=list()
for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}
#绘制箱线图
#set color

ggplot(data = data,aes(x = clinical, #分组列名
                       y = gene, #连续变量列名
                       fill = clinical))+ #按分组填充颜色
  scale_fill_manual(values = ggsci::pal_npg("nrc", alpha = 0.6)(9)) + #用自定义颜色填充
  geom_violin(alpha = 0.4, position = position_dodge(width = .75), 
              size = 0.8, color="black") +
  geom_boxplot(notch = TRUE, outlier.size = -1, 
               color="black", lwd=0.8, alpha = 0.7) +
  geom_point(shape = 21, size=2, # 点的性状和大小
             position = position_jitterdodge(), # 让点散开
             color="black", alpha = 1) +
  theme_classic() + 
  ylab(gene) +
  xlab("") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
        #panel.border = element_rect(colour = "black", fill=NA, size=0.2),
        axis.ticks = element_line(size=0.2, color="black"),
        axis.ticks.length = unit(0.2, "cm"),
        legend.position = "none",
        axis.title = element_text(size = 15),
        axis.text = element_text(size = 10)) +
  theme(panel.border = element_rect(fill=NA,
                                    color="black", 
                                    size=1, 
                                    linetype="solid"))+
  # 如果不要组间比较就注释掉下面这行
  stat_compare_means(comparisons = my_comparisons,
                     label.y = max(data$gene)-0.05)
ggsave(paste0(gene,".pdf"), width = 4, height = 4)

